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We extend the Collins-Soper-Sterman (CSS) formalism to apply it to the spin-dependence gov- 
erned by the Sivers function. We use it to give a correct numerical QCD evolution of existing 
fixed-scale fits of the Sivers function. With the aid of approximations useful for the nonperturba- 
tive region, we present the results as parametrizations of a Gaussian form in transverse-momentum 
space, rather than in the Fourier conjugate transverse coordinate space normally used in the CSS 
formalism. They are specifically valid at small transverse momentum. Since evolution has been 
applied, our results can be used to make predictions for Drell-Yan and semi-inclusive deep inelastic 
scattering at energies different from those where the original fits were made. Our evolved functions 
are of a form that they can be used in the same parton-model factorization formulas as used in the 
original fits, but now with a predicted scale dependence in the fit parameters. We also present a 
method by which our evolved functions can be corrected to allow for twist-3 contributions at large 
parton transverse momentum. 



I. INTRODUCTION 

High energy collisions with transversely polarized 
hadrons are ideal processes for extracting information 
about the structure of hadrons. The nonperturbative 
functions that enter into the corresponding factorization 
formulas are sensitive to novel aspects of QCD dynamics 
such as chiral symmetry breaking and the role of orbital 
angular momentum. (See e.g. [l| for some interesting 
recent discussions.) The Sivers function is an example 
which has received considerable attention in recent years, 
and will be the focus of this article, although many of 
the results and techniques are extendable to other in- 
teresting transverse-momentum dependent (TMD) func- 
tions. In loose terms, the Sivers function describes the 
transverse-momentum distribution of (unpolarized) par- 
tons inside a transversely polarized hadron (usually a 
proton). In semi- inclusive cross sections with a single 
transversely polarized target hadron, it leads to a char- 
acteristic sin((/) — azimuthal modulation (0 and (jfh 
being the azimuthal angles of the transverse spin and the 
produced hadron, respectively). It is one of a collection 
of TMD parton distribution functions (PDFs) and frag- 
mentation functions (FFs) that are actively being studied 
for the insight they can provide about hadron structure 
and the unique opportunities they provide for comparing 
theoretical descriptions with experimental results [2t6|. 

The Sivers effect was originally proposed more than 
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two decades ago in Ref. as a mechanism for gen- 
erating transverse single spin asymmetries (SSAs) in 
hadron-hadron collisions. Shortly afterward, it was ar- 
gued in Ref. 8] on the basis of time-reversal (actually 
TP) invariance that the Sivers function vanishes. This 
result, if true, implies that the corresponding SSA in 
semi-inclusive deep inelastic scattering (SIDIS) is power- 
suppressed (i.e., it is of "higher twist"), leaving only the 
spin-dependent effects due to the Collins function in frag- 
mentation. Thus, a contradiction arose when specta- 
tor model calculations [9] gave an explicit nonvanishing 
leading-twist SSA in SIDIS with the azimuthal depen- 
dence associated with the Sivers function. The situa- 
tion was clarified in Ref. where it was shown that 
the proof of vanishing of the Sivers function was incor- 
rect in QCD, because it ignored the Wilson lines needed 
in the definitions of parton densities. Instead, the true 
consequence of TP invariance of QCD is that the Sivers 
function reverses sign between SIDIS and Drell-Yan (DY) 
processes. This is because future-pointing Wilson lines 
are needed in TMD functions like the Sivers function 
when used for SIDIS, but past-pointing Wilson lines are 
needed for the Drell-Yan process. At the level of the 
actual cross section, the sign-reversal for the Drell-Yan 
process was verified in model calculations in Ref. (l]| . 

Certain other polarization or azimuthally dependent 
functions, such as the Boer-Mulders and the pretzelosity 
distributions [H, US]: also share this "T-odd" property 
of reversal of sign between SIDIS and Drell-Yan. Over 
the past decade, there has developed much work in the 
extraction, study, and formal theoretical description of 
these functions. 

However, phenomenological fits of the Sivers function 
(and of related functions) have so far [3, [l^ used only 
the simplest parton-model factorization formulas where 
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the TMD parton densities and fragmentation functions 
do not evolve with the scale of the process, or use incor- 
rect evolution formalisms. This is inadequate when they 
are to be applied to experiments at widely different ener- 
gies. There is a good QCD formalism for applying TMD 
functions in a factorization framework, due to Collins, 
Soper and Sterman (CSS) [11,113. The CSS formahsm 
gives a correct treatment of the region of low transverse 
momentum, which is where the Sivers function analysis 
is used. However it has not been fully systematized for 
the case of the Sivers function and other azimuthally- 
dependent functions, except in the work of Boer [l^ [l^ 
and Idilbi et al. ^] , on which we comment below. 

In this paper, we give a complete extension of the CSS 
method to processes that need the Sivers function, using 
the methods recently given in Ref. It is straight- 

forward to extend our results to the other azimuthally 
dependent PDFs and FFs (e.g., the Collins function and 
the Boer-Mulders function). We apply the formalism to 
give numerical results for the Sivers function evolved from 
existing fits. The only extra nonperturbative information 
needed for the evolution is universal and is obtained from 
existing fits to the unpolarized Drell-Yan process. This 
extends the results given by two of us in Ref. [l^ for 
the unpolarized case. Reference [l^ attempts to include 
some effects of evolution by simply including the evolu- 
tion from collinear factorization, but this is incorrect for 
TMD-factorization. It is also stated (Ref. [l^ , for exam- 
ple) that the true scale-evolution of the Sivers function is 
unknown. One purpose of this article is to demonstrate 
that this is no longer true. 

With the aid of an approximation useful for the non- 
perturbative region, we present the results as Gaus- 
sian transverse-momentum distributions with scale- 
dependent parameters. They are therefore as easy to use 
in simple parton-model-style calculations as the original 
fixed-scale fits [13, EBl- As the scale increases, the distri- 
butions broaden substantially in transverse momentum, 
and get diluted in size. It will be necessary to include 
perturbative twist-3 corrections to get more accurate val- 
ues at the larger values of transverse momentum, and we 
present a scheme for how this should be done. 

Boer [l^ [3 has applied the CSS method to processes 
involving the Collins function. Idilbi et al. [l^l have ap- 
plied the CSS method to their definitions of various TMD 
distributions [H, including the Sivers function. Our 
treatment is substantially improved, to include a correct 
treatment of the nonperturbative region in CSS evolu- 
tion applied to T-odd functions, to use a more modern 



version of the CSS formalism, to apply it to the Sivers 
function, and to obtain convenient numerical results for 
the Sivers function. 

Although it has recently become common for the 
word "resummation" to be used to indicate any CSS- 
like treatment, in our work we will maintain a firm dis- 
tinction between resummation methodology and TMD- 
factorization. The term "resummation" is often used to 
indicate that one starts with conventional collinear fac- 
torization and resums logarithms of qr/Q, which can in 
fact be done with the CSS methodology. The problem 
with this approach is that it is only valid when the un- 
derlying collinear factorization formula is valid, i.e., for 
the region where the transverse momentum qt is both 
much less than the hard scale and much greater than 
hadronic binding energies ~ Aqcd- (See, in particular, 
the recent work of Ref. [Hf.) But to extend the cal- 
culations to transverse momenta comparable to Aqcd 
and to zero transverse momentum requires a complete 
TMD-factorization formalism, which we use here. This is 
particularly important because many SIDIS experiments 
such as HERMES and JLab are performed at kinemati- 
cal scales where transverse momenta of order Aqcd are 
certainly important, and Q is not so large. 

A number of difficulties are caused by the use of a pure 
resummation formalism rather than TMD factorization 
as the basis of calculations. For the present paper, one of 
the most significant is that a leading-power resummation 
formalism does not give the effects associated with the 
Sivers function (and also those associated with the Boer- 
Mulders [111 function). But, provided that spin effects 
are treated correctly, the presence of these functions is 
automatic in TMD factorization, at leading power. 



II. SETUP AND DEFINITIONS 

In this section we give the factorization formula for 
SIDIS: e + P{S) — > e + h + X, and present the definitions 
of the TMD functions. We let P and S be the momentum 
and spin vector of the hadron target, and we let h label 
the detected hadron, of momentum ph ■ With a single ex- 
changed photon of momentum q, independent kinematic 
variables are: Q — \/— g^, x = Q'^/2p-q, z = P-ph/P-q, 
and the virtual photon's transverse momentum q-p (in 
a hadron frame where the measured hadrons have zero 
transverse momentum). 



I 

The TMD-factorization formula in the form derived by Collins [2l| is: 

+ r(Q,qT) + o((A/g)'^). (1) 
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Here F^/pt (a;, kiT, 5") is the TMD PDF for an unpolar- 
ized quark of flavor / in a proton of polarization S, and 
Dhif{z, zk2T) is the unpolarized fragmentation function. 
These factors contain nonkinematic parameters, /i, Cf, 
and C-D: whose definitions are given below. The hard- 
scattering factor I'H^I'"' is computed, with appropriate 
subtractions, from massless parton scattering in a photon 
frame where the photon and partons have zero transverse 
momentum — see [111, page 527] for its definition. The 
first line of the factorization formula is valid at low trans- 
verse momentum, and the Y term provides a correction 
for large transverse momentum in a form like that for or- 
dinary collinear factorization. Although we will focus on 
SIDIS for this paper, the same general treatment applies 
also to DY scattering, up to the change in direction of the 
Wilson line in the definition of the TMD PDF. Note that 
the TMD- factorization piece, the first term in Eq. ([IJ, is 
formulated specifically to deal with the small fc^ behav- 
ior (fcy — 7> 0), while allowing for systematic corrections 
to the behavior as kr grows larger than Aqcd- 

The above formula is exactly like the parton-model 
formula for the same cross section except for the scale 
dependence of the PDF and fragmentation function and 
except for higher-order corrections in the hard scattering 
and y-term. It differs from the older CSS formula by no 
longer needing an explicit soft factor. The factorization 
formula ([IJ is written for the case that the partons at 
the hard scattering are unpolarized. Parton polarization 
effects can be allowed for simply by inserting spin ma- 
trices for the incoming and outgoing partons of the hard 
scattering. This gives other terms, e.g., those with the 
Collins function in fragmentation, with their character- 
istic angular distributions in the cross section. It was 
recently suggested in Ref. [J?! that it would be useful 
to analyze data for cross sections in transverse coordi- 
nate space hr by taking various weighted integrals with 
Bessel functions. In that case, the &t version of Eq. ^ 
is needed. 

The parameter is a conventional renormalization 
scale, which we will choose to be in the MS scheme. 
It should be chosen to be of order Q so that the hard 
scattering has no large logarithms. The parameters C-F 
and C,D are related to the need to regulate rapidity di- 
vergences in the definitions of the TMDs. They are de- 
fined with the aid of an auxiliary rapidity parameter , 
which has the function of separating forward and back- 
ward rapidity gluons. We use a hadron frame (in which 
the hadrons have zero transverse momentum), oriented 
so that e^^ ^ e^^h , and we let Mp and Mh be the masses 
of these hadrons. Then C^p and Cd are defined by 



and 



(3) 



They obey \/CfQd — up to power-suppressed correc- 
tions, and have been normalized to correspond to CSS's 
definitions. 

The definitions of gauge-invariant TMD functions are 
equipped with Wilson lines. A Wilson line (or gauge link) 
from a point a: to cx) along the direction of a four- vector 
n is defined as 



W(po, x\n) = P exp 



poo 

I ds n ■ Al{x + snY' 
Jo 



(4) 

Here, bare field operators and bare couplings are used 
and P is a path-ordering operation. The generator for 
the gauge group in the fundamental representation, with 
color index a, is denoted by t"' . 

To define the parton densities, we use two lightlike di- 
rections that characterize the extreme forward and back- 
ward directions: 



ma = (1,0,Ot), mb = (0,1,Ot). 



(5) 



These correspond to the directions of P and ph- Our 
coordinates for a 4-vector v are defined by 



V = {v'^ , V , vt) 



where. 



(6) 



(7) 



Now the most obvious definitions of PDFs use light- 
like Wilson lines, which give rise to rapidity divergences 
f^S^]. Regulating the divergences can be done by using 
non-light-like Wilson lines. So we define vectors riAii/A) 
and riB(2/B) with finite rapidities jja and ys'- 



riA = (1, 



riB 



(-e2^^1,0T). (8) 



2„2 2{vp-y,) 



(2) 



The actual TMD PDF in Eq. (P is defined as a hmit 
of an unsubtracted TMD multiplied by certain unsub- 
tracted soft functions. These are first defined in trans- 
verse coordinate space and then the final result will be 
Fourier transformed to transverse-momentum space. The 
unsubtracted TMD PDF is 
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F;-f(a;,bT,5;A^;yp-2/i3) 

^TrcTrz, J ^e-"''^'^' {P, S\i^f{w/2)W{w/2,oo,nB{yB) f^W{~w/2,oo,nB{yB))i^fi~w/2)^^^ (9) 

where w'^ = (0^, , hx), and we notate the functions with a tilde to indicate the use of transverse coordinate space. 
The subscript c indicates that only connected diagrams are included, and Trc and Tr/^ represent color and Dirac 
traces respectively. The unsubtracted soft function is 

5(o)(bT;yA,yB) = ^(0|VK(bT/2,oo;nB)L VK(bT/2,oo;nA)adW^(-bT/2,oo;nB)bcW^(-bT/2,oo;nA)L|0). (10) 

In both of these functions, there should be inserted transverse gauge links at infinity. However, their effects cancel in 
the subtracted TMD PDF, when Feynman gauge is used, so we have not indicated the extra gauge links explicitly. 
The fuU definition of the TMD PDF from [21] is 



Ff^p,{x,hT,S;^i,CF) = FJfprix,hT,S;^i■,yp^i~^)) ^(o)(bT; +oo, y,) ^^^^ 

^ i(o)(bT; +0O, -cx))6(o) (bx; 2/s, -oo) 

This involves limits: infinite rapidity on the Wilson lines indicated, infinite length for the Wilson lines, and then 
removal of the UV regulator (dimensional regularization). The factors ZpZ-i at the end of Eq. ([TT]) are the field 
strength and TMD renormalization factors respectively. Notice that two of the soft factors have one of their rapidity 
arguments equal to the finite parameter y^. 

An exactly analogous definition applies to the fragmentation function (see Ref. [2lj for the explicit definition). In 
our notation, capital letters will denote unintegrated quantities and lower case letters will denote quantities integrated 
over transverse momentum. Otherwise, we will stick as closely as possible to the Trento conventions p9| . 

The momentum-space TMD PDF is 

F//Pt(x,kT,5;/i,CF) = ^ ^ ShT^'^^'"^ F^IP,(x,h^,S-iiXF). (12) 

This has dependence on the azimuthal angle between and the transverse spin vector St of the target hadron. (We 
normalize St so that its maximum size is unity.) The TMD PDF is decomposed as usual into the unpolarized TMD 
PDF and a spin-dependent term: 

Ff/pt{x,kT,S;^i,CF) = Ff|p{x,kT]^i,CF) - F^j.^ {x, kp; ^J., Cf) ^'^'^^ , (13) 

with F:^rjl (x, kp] fi, (^p) being the Sivers function. 



III. EVOLUTION OF THE SIVERS FUNCTION A. Coordinate Space Representation of Azimuthal 

Dependence 



In this section we generalize CSS evolution from the 
unpolarized TMDs to the Sivers function. Similar meth- 
ods apply to the other TMDs with azimuthal dependence. 

The general CSS formalism works equally well for these 
functions [21]. But it involves Fourier transformations in 
two transverse dimensions, and for practical use it is con- 
venient to perform the azimuthal integrals analytically 
and to write the transforms in terms of integrals over the 
sizes of the transverse variables. The treatment of the az- 
imuthal integrals provided in Sec. | III Al closelv parallels 
previous treatments in Refs. [20|.|23| a"iid recently in [27| . 



To analyze the evolution of the last term in Eq. (fT3)) 
we extract the azimuth-dependent part by defining 

0}/p(x,kT;Ai,CF) = -^F^J{x,kp-^iXF), (14) 
in terms of which the complete Sivers term is 

F^r^^{x,kp;fi,CF) '•^^J ^ = 0}/p(x,kT;Ai,CF)eyS'T- 

(15) 
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The Fourier transform of the Sivers function is 
F^/{x, hr; Cf) = j d'kT e"^''-'^- F^J {x, kT;fi, Cf) 

= 27r / dkTkTJo{kTbT)F:^T{x,kT;fi,CF), (16) 
Jo 

and the Fourier transform of (l)^j^p{x, kx; /x, C,f) is 



p(x,bT;/i,CF) - / d2kTe-''^^-^^0}/p(x, kT ;m,Cf) 



d kx e 



F^jJ [x^kT^nXp) 



(17) 



Using Eq. (|T6|) gives 



1 6* ~ 

(a;,bT;Ai,CF) = * ■P'i't '^(a;, &t; Ai, Cf) , (18) 



where we have denoted the derivative of F^j,-' with re- 



spect to the length of b-r by 

_ dF^f{x,bT;fJ.XF) 



F[^^{x,bT;n,CF) 



dh 



(19) 



As we will see shortly, it is this derivative F and not 
the function F itself that gets used in the evolution equa- 
tions and in the formula for the Sivers term in the actual 
transverse- momentum dependence in Eq. (IT^ . 

Taking an inverse Fourier transform of Eq. (|18l) allows 
(f)^j^p{x, kx; fJ,, Cf) to be rewritten in terms of Eq. (fK 

(?!)}/p(a;,kT;M,CF) 

^ J d'hT e*-''- ^}/p(.T,bT;A^,CF) 



(27 



(27r)2Afi: 



d^bT e'''-'^-p^i?i'^/(a;,6T;M,CF) 



(20) 



To further simplify this expression, and without loss of 
generality, we use a frame where kx is in the x direction 

so that = (1,0) and = (cos0,sin0). Then, 



'f/P 



(a;,kx;M,CF) 



(2^)2Afp Jo 
(2^)2Mp / 



dbTbTF[^f{x,bT]^iXF) / d6le"^^''^™''''(cos6l,sin 



Ak-T^T COS d 



dbx br F[;^^ [x, br; fi, Cf) 



d 



d{kTbj 
d 



^g, g.feTbT cos «(-^^Q^ 



■MkTbT) 



2TTMpkT Jo 

Then the complete Sivers term in Eq. (|13l) is 



d&T hTJi{kTbT)F[^^ [x, br] fi, Cf) 



(/)}/p(a;,kx;Ai,CF)eij5'^ 



^ ^ ' dbTbTJlikTbT)F[^^ix,bT;tlXF) 



271 MpkT Jo 

So, from Eq. (jlSp we express the momentum-space Sivers function in terms of F 



F^r^^{x,kT;fi,CF) = 



-1 
27rfcT Jo 



dbr bTJi{kTbT)F[^^{x, br; n, Cf) 



whose inverse transform is 



-F'lT ■''(a;, &t;M, Cf) = -27r / dkr k^Ji{kTbT)F^j/ {x,kT; fi,CF 



(21) 



(22) 



(23) 



(24) 



Notice that the originally defined F^-^ from Eq. (|T6l) no longer appears. The ^T-dependent function F[^^ {x, &t; Cf) 
is closely analogous to the quantity /^^^^ that appears in Eqs. (16) and (20) of Ref. [13], and to dlqr in Eq. (40) of 
Ref. [20], though the basic definition for the fo^-space TMD PDF in Eq. ([TT|) is significantly different. 



B. The Evolution Equations 



to Cf, and the RG equations which give evolution with 



The set of evolution equations comprises the Collins- 
Soper (CS) equation which gives evolution with respect 



6 



respect to fi. The CS equation for the TMD function 
defined in Eq. (HH) is [U 

dFf/pt{x,hT,S;n,CF) _ 



KibT;ti)Ff^prix,hT,S;^iXF), (25) 



where 



f^fu \ 1 ^ 1 I S{bT;ys, -oo) 



2 \S{hT]+(yo,ys) ^ 

The RG equations are 



din /I 

and 

dFf/p^{x,hT,S]^i,CF) 



'lK{g{p)) 



(26) 



(27) 



din 



= 7i=^(.9(M); CF^')^^//Pt (a:, bx, 5; A^, Ci=^). (28) 

Similar equations apply to the fragmentation function. 
It follows that the Qp dependence of 7f is determined: 



9 In vCf 



= -iKigi^J-)), 



(29) 



SO that 



1 Cf 
7F(ff(Ai);CF//i^) = 7F(g(A*);l) - 27-R'(3(A*))ln^. (30) 

These equations were used in Ref. [l^] to calculate 
the evolution of the unpolarized TMDs. For the spin- 
dependent case, the Fourier transform of the second term 
in Eq. ([T3l) obeys the same evolution equations, i.e., the 
equations apply to 



J ci2kTe-^''-b-i^^iV(x,fcT;M,CF) 



= 4>)/p{x,hT■,^J.,CF)e^JS^p. (31) 
The CS equation for the spin-dependent part is therefore 

c>i>)/p{x, bx; A^, C,F)tijSip 



91n vCf 

= K{bT; n)(t>}/p{x,hT; fJ., CF)£ijS^j, . (32) 



Hence, Eq. ([T8| shows that the CS equation for 
^ {x,bT', fJ:, Cf) is the same as for the unpolarized 



TMD PDF: 



d\nF[^^{x,bT;fi,CF) 
91n VCf 



KibT;^i). (33) 



Similarly, its RG equation is like Eq. ([28| : 



dF^T^{x,bT;n,CF) 



dlnfi 



= 7f(5(m); CF/^i^)F[^Hx, bp] ^l, Cf). (34) 



Note that in Eqs. ^3} and ((Ml) the same CS kernel 
K(bT]lJi) and anomalous dimension jF{g{p)', Cf / IJ^'^) ap- 
pear as in the unpolarized case. This is because K and 'jp 
are properties of the operator defining the parton density, 
and this operator is the same for the ordinary unpolar- 
ized TMD PDF as for the Sivers function; both concern 
the number density of quarks in a hadron, with no po- 
larization restriction on the quark. 

It is important to emphasize that the evolution equa- 
tions (PSI \T7\ [25)) are set up to be exactly correct for 
all bp^ and for all kp. This includes the region where 
6t — oo (and hence kp 0). Indeed, the first term on 
the right side of Eq. ^ (the TMD-factorization term) 
is designed to give an accurate pQCD treatment when 
kp Q, independently of the relative sizes of kp and 
Aqcd- 

C. Power laws for kp and bp dependence 

As a guide to the qualitative behavior of the Sivers 
function, we summarize in this section the power laws 
for its dependence on transverse momentum and trans- 
verse position as obtained from simple model calcula- 
tions. (For a detailed treatment of the power law behav- 
ior of other TMDs, see Ref. (S^l and also recent discus- 
sions in Ref. (27j.) In purely perturbative higher-order 
calculations, these get modified by logarithms, while use 
of a correct solution of the evolution equations can sig- 
nificantly modify the power laws [3lj. Nevertheless, the 
power laws from elementary perturbative calculations 
form a useful standard of comparison. 

First, we characterize the power law for an ordinary 
unpolarized TMD PDF by 



F{x,kq 



1 



(35) 



At large kp, the falloff 1/kp is the simple dimensional- 
analysis power, appropriate to a theory with a dimen- 
sionless coupling. The increase at low kp is tamed by an 
infra-red cutofi^ M, which in QCD is nonperturbative. In 
bp space, the large-kp behavior Fourier transforms to 

F{x, bp) ^ constant x logarithms (as bp 0). (36) 

At large bp, the falloff of F should be at least rapid 
enough that the integral over all bp is convergent, to give 
a finite value for F{x,kp) at kp = 0. Normally an expo- 
nential or Gaussian falloff is assumed (which is controlled 
by nonperturbative effects in QCD). 

As for the Sivers function, its contribution to the quark 
density , F^p^ {x, kp)eijkpS^ / Mp, has a kinematic zero at 
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kx = 0. In addition, it is a chirality- violating quantity, 
and at large fcr j this requires a suppression by a factor of 
mass divided by kr relative to the unpolarized density. 
So we characterize the result by 



krM 



(4 + M2)2- 



For the Sivers function itself, we therefore have 

,^f, . , M2 



F^r/{x, kr ) 



+ M2)2 



(37) 



(38) 



This falloffis characterized as "twist-3." In hx space, the 
behavior of the Fourier transform of (p8l) at small 6t is 



X, bx) ^ constant + 5^ x logarithms. 



(39) 



However, as we saw, it is the derivative of this quantity 
with respect to hx that is actually used, for which the 
behavior is linear: 



6t X logarithms. 



(40) 



Although the actual equations for evolution are the same 
for the Sivers function as for the standard unpolarized 
TMD PDF, there are substantial differences in the way in 
which the evolution is reflected in the numerical values of 
these functions in transverse-momentum space. Because 
F-^^ is approximately linear in bx at small 6t and because 
the Ji Bessel function instead of Jq appears in Eq. (21), 
the Fourier transform for the Sivers function is sensitive 
to larger bx values than the transform for the unpolarized 
TMD. This also implies that the evolution of Sivers is 
subject to more uncertainty from the nonperturbative 
large-^T region than that of the unpolarized TMD. 



D. Small-^T expansion 

For the unpolarized TMD PDF, an expansion for small 
br can be made in terms of the integrated PDFs. After 
Fourier transformation, this gives both the large-fc^ be- 
havior, and the normalization of the integral over the 
whole small kx region. 

The same idea continues to apply when we include the 
dependence of the TMD density on the target polariza- 
tion. We can write 



F(x, hT,S) — ^ coefScientj 



(P,S'|operator|P,S'), 



(41) 

where the coefficients and operators are unaltered since 
they are properties of the TMD number-density operator. 
But the twist-2 operator on the right-hand side of (|41l) is 
the ordinary number-density operator used to define an 
integrated PDF, and its matrix element is independent of 
transverse spin. Thus the twist-2 operator, correspond- 
ing to a l/fc|n fall off at large /ct, provides no contribution 
to the Sivers function in Eq. (|411) . The leading large-fc^ 
behavior of the Sivers function is the 1/k^ term associ- 
ated with the twist-3 operators, the same operators that 
are used in the Qiu-Sterman formalism i 32|] . 



IV. OBTAINING EVOLVED SIVERS 
FUNCTIONS 

In this section, we discuss the steps for obtaining the 
evolved Sivers function using already existing fits to the 
nonperturbative parts. 



A. Solution in terms of fixed-scale Sivers function 

Previous fits p^ . [Tsj of the Sivers function used the 
parton-model formula for the hadronic tensor. We now 
show how these can be converted to use the correct QCD 
formula. 

The parton-model version of TMD factorization 
amounts to applying the following approximations to the 
true QCD formula (H)): 

(i) Replace the hard scattering by its lowest order. 

(ii) Neglect the Y term. 

(in) Omit the evolution of the TMD PDFs. 

If the renormalization scale /i is taken of order Q, higher- 
order corrections to the hard scattering are purely pertur- 
bative. One of the simplifications for TMD factorization 
is that these are just an overall factor, dependent on Q 
only through the running coupling as{Q)- This factor 
is the same, independent of the hadron and the quark 
polarization, so it does not affect the ratio of the Sivers 
function to the ordinary TMD PDF. 

The Y term only affects large transverse momentum 
(of order Q), whereas the data is dominantly at trans- 
verse momenta in the nonperturbative region. So the 
neglect of Y should be an adequate approximation with 
present data, and is easily corrected in the future, with 
the aid of fits for the Qiu-Sterman twist-3 function. 

For a fixed value of Q, the TMD functions can be given 
fixed values of /i and C,f, ^i = Q and Cf — Q^, and the 
QCD factorization formula is the same as the parton- 
model formula, up to an overall if-factor. This legit- 
imizes the fixed-scale fits. But as can be seen from Fig. 
[T] below, evolution gives substantial changes in the TMD 
PDFs needed at higher Q. These are easily obtained, in 
their transverse-coordinate-space form, in terms of the 
parton-model fits at a fixed scale. We derive the neces- 
sary result starting from Eqs. ([55)1 . ([Ml) , and ([501) . 

In these equations, the anomalous dimensions and 
7k are perturbatively calculable, but the function K at 
large values of bx is nonperturbative. We follow Ref. [l3| 
to separate the perturbative and nonperturbative parts 
of K. First, we define 



fJ-b 



0* 



Here Ci is a fixed numerical coefficient and timax is 
chosen to keep in the perturbative region. In the 
fits to unpolarized Drell-Yan, the values chosen were 
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b^^^ = 0.5 GeV-i in [H, and W = 1-5 GeV"! in 
Next we write 

K{bT;fi)^K{b,;fib)~ r %iMli'))-9K{hT). (43) 

The first two terms are perturbative and include all the 
evolution of K. The last term is nonperturbative but 



scale independent. It represents the only nonperturba- 
tive information needed to evolve the Sivers function from 
the scale Qo where it was initially fit. But this function is 
process independent f^Tl , so we can take its value from al- 
ready existing fits to unpolarized Drell-Yan (33i . i34i | scat- 
tering at a variety of energies. 

This gives the evolved function: 



' 



7F(g(A*'); 1) - In — r^K{g{^^')) 



— In -p^lKigifJ' )) - 9K{bT) In — — 
^0 M Qo Qo 



(44) 



We can set fiQ = Qq and then use Qo = V2.4 GeV, which 
is the appropriate scale for the fits in [3: 113 , which used 
data from the HERMES experiment. For the prediction 
of data at a higher energy, one should set ji^ = C,f = 

. The anomalous dimensions 7^ and 7^ are used in a 
region where perturbative calculations are appropriate. 

The Sivers function in transverse-momentum space is 
then obtained from Eq. (|44l) by Fourier transformation, 
as in Eq. ([23l) . 

The one-loop values of the relevant perturbative quan- 
tities are listed in the Appendix. 

The size of the Sivers asymmetry is also often 
parametrized by the function 



Ff|p^(x,\^.^;\S,\xXF) - -F^z/pt (a;, kx; -5, ^, Cf) 
where 

A^%pt(x,fcT) = -^-^F^J^xM-.P^Xf). (46) 



As can be seen from Figs. [T] and [2] below, TMD func- 
tions broaden substantially as the scale increases. Thus 
larger values of transverse momentum become important, 
and correspondingly we need the F factor at small 6t- 



Including the perturbative calculation of Sivers 
function at small-br 



At low scales, the Sivers function is dominantly at low 
values of fc^, and correspondingly the range of 6t that 
matters concerns the larger values where both the start- 
ing value ^(x, 6t; MOj Qo) *nd the evolution kernel 
K{bT', ^J■) are in the nonperturbative region. After evolu- 
tion to a sufficiently large scale, the broadening of the fc-r 
distribution makes smaller values of bx important, where 
there is perturbative information. For both this case and 
the treatment of the large-Zc^ tail of the Sivers function 
we can use the expansion (|4ip to write it in terms of the 
twist-3 Qiu-Stcrman function. 



Following the method used for the unpolarized TMD PDF — see Ref. [T3,[2H and Eq. (31) of Ref. [221 — we write 

Fl^^^{x,bT;fJ.XF) = X! ^^2^^ / CjY°''%xi,X2 , b*;^il,^ib,g{fJ-b)) Tp j/p{xi , X2, fJ-b) 



X Xi Xi 



xexp(ln^7^(6,;M.)+ T ^' 



iF{g{p^); 1) - In — —iK{g{p^)) 



I xexp|-5Siv-(^,fe^)_g^(;,^) 



In 



Qo 



(47) 



The first line describes the matching to a 
coUinear treatment relevant to small b^. There, 
F-^^^ {x,bT', (J-tCf) is expressed as a coefficient function 



C f /j{xi,X2,b^:; ^1, ^,b,g{^b)) convoluted with a (twist-3) 
Qiu-Sterman function Tp j^p(xi, £2, fJ-b), where for the 
simplicity, we neglected the terms proportional to the 
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derivative of the twist-3 Qiu-Sterman function. On 
the second hne, the first exponential comes from the 
perturbative part of the evolution of the Sivers function; 
the use of 6* and fib ensures that C, K, jp, and "/k are 
in the perturbative region. The second exponential gives 
a correction to allow for nonperturbative behavior at 
larger bx- In its exponent are both the nonperturbative 
term gxibr) for the evolution kernel, and an extra term 
g^Jp^ix, bp) for the Sivers function itself. These terms 
are both scale independent. ^ 

The coefficient C can be determined, for exam- 
ple, by performing a low-order perturbative calcula- 
tion of the left-hand side of Eq. (|T7)) . of the Qiu- 
Sterman function, and of the first exponential, while 
ignoring the nonperturbative correction [251 . The nor- 
malization factor, Mpbp/2, in Eq. (|T7)) ensures that 
Tp{xi,X2, fib) has the standard normalization [2^, and 
at zeroth order the contribution to the hard coefficient 



IS 



Sivers, (0) 

f/i 



{xi,X2,b^^.;nl,Hb,g{pb)) 



Sf^jSil-x/xi)S{l-x/x2), (48) 



which is similar to the zeroth order term in Eq. (All) 
of Ref. [12] for the unpolarized case. (Recall that, since 
the Qiu-Sterman function Tp{xi,X2, fJ-b) is universal, an 
extra minus sign is needed if we consider Drell-Yan in- 
stead of SIDIS.) The factor of bp in the normalization 
is a reminder that it is the derivative of the Sivers func- 
tion that we evolve in Eq. (H7|) . not the Sivers function 
itself. Higher-order contributions to the coefficient func- 
tion can be taken directly from work, such as Ref. [25| . 
which treats smaller bp within the Qiu-Sterman method. 
Calculations of the unpolarized coefficient functions to 
higher orders in the MS scheme have already been car- 
ried out in Ref. 

The corresponding formula for the unpolarized TMD 
PDFs is very useful, since instead of the Qiu-Sterman 
function it uses the ordinary integrated PDFs, which are 
very well measured. In contrast, the phenomenology of 
the Qiu-Sterman function is less well known quantita- 
tively, so there may be less of an advantage of using Eq. 
(|47l) instead of Eq. (|44| . 

In the remaining sections, we will discuss the imple- 
mentation of evolution, given some nonperturbative in- 
put functions, and provide specific evolved fits. Before 
continuing, however, we should emphasize that matters 
related to the fitting of the nonperturbative functions, 
including the choice of functional form for gK{b) and 
the matching procedure in Eq. (|42|) . are unrelated to 
the validity of the TMD-factorization formalism itself. 
The TMD-factorization formalism automatically accom- 



Note that our sign convention on g^Yp^(x,bT) and gKibp) is 
opposite of Ref. [5^ . 



modates any refinements to knowledge about the non- 
perturbative physics. Indeed, a central aim of this arti- 
cle is to demonstrate the generality of the method. In 
our calculations below, we have chosen to consider fits 
to the nonperturbative functions that correspond to de- 
tailed studies of existing data. In addition to providing 
tools for phenomenology, our calculations illustrate how 
numerical values for the Sivers function corresponding to 
the definition in Eq. (jlip can be obtained, once the non- 
perturbative functions are constrained by data. Thus, 
our use of TMD-factorization is closely analogous to what 
already exists for collinear factorization. 



V. GAUSSIAN PARAMETRIZATIONS IN THE 
LOW-gT REGION 

In this section we explain the implementation of QCD 
evolution for the Sivers function with a Gaussian ansatz. 
Since the small-^T region is twist-3, the tail of the 
(momentum-space) Sivers function (at large kp) is power 
suppressed relative to the unpolarized TMD function. 
Furthermore, as illustrated in Ref. [22], a Gaussian 
parametrization provides a good description of the low 
transverse-momentum behavior, even up to transverse 
momenta of a few GeV. Therefore, we take as a start- 
ing point a detailed treatment of the twist-2 large-fey be- 
havior, leaving for future refinements an account of the 
matching of the small-bp behavior to the twist-3 factor- 
ization formalism. That is, we use Eq. (j44p rather than 
Eq. gZl) 

Even so, a full treatment that extends to small-fer by 
including higher orders in Cf /j{xi,X2,bf,; fj-l, fj,b, g{fib)) 
will be crucial in the long run for a complete understand- 
ing of the evolved Sivers function over the full range of bp . 
This is especially important to keep in mind when deal- 
ing with weighted integrals of the Sivers function where 
the effect of the large transverse-momentum tail becomes 
magnified. We intend to pursue this in future refinements 
of the TMD approach. 

At the initial fitting scale, we drop the explicit scale 
dependence: 

Fi^^oi^^bp) = F^^f {x,bp; fio,Ql)- (49) 

To match previous fits [l^.lT5|. we approximate the input 
function by a Gaussian 



F[^Ux,bp) 



{kp)aftp{x)b'i 



exp [^{kDob^p/A] , 

(50) 

which corresponds also to a Gaussian ansatz for the 
momentum-space distribution: 



Fip^Qf{x,kp) = 



/iV(^ ) 



exp 



T/Q 



(51) 



The parameter (fc|i)Qis the width of the Sivers function 
for a quark of fiavor / at the scale where the Gaus- 
sian fit is performed. Comparing with Eq. (j47p . we 
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see that gfj'^'%x,bT) = (fc2)^6|/4. The fits performed 
in 1^ [H] are for quite low scales (Q^ = 2.4 GeV^ for 
HERMES data). We therefore assume that the Sivers 
function is dominated by the nonperturbative large-&T 
region, in which case a Gaussian description, with a neg- 
ligible tail effect, makes sense. The first moment of the 
input momentum-space Sivers function obeys the usual 
relation; 



/: 



it! 0^ (''') 



2M, 



(52) 

We again remind the reader that for our calculations, we 
are assuming a Sivers function for SIDIS and that a sign 
flip is necessary to go to DY. 

With gKibx) already known from previous fits to high 
energy Drell-Yan data [ssl - lssj . all that is now needed 
in order to obtain evolved Gaussian fits are (fc|-)o and 
/j^^ (x) . These will come from previously obtained fixed- 
scale Gaussian fits. In the next section, we will provide 
two examples and illustrate the effect of evolution for two 
of the sets of Gaussian fits available in the literature. 

The function 5_r-(6t) is the only nonperturbative input 
that is necessary apart from these initial fits. We have 
also adopted the standard Gaussian ansatz for gjf (6t) , 
writing gKibr) — 52&t/2. Fits like those of Refs. |33l - [35| 
provide numerical values for 52- In the Brock-Landry- 
Nadolsky-Yuan fits [H a value of 52 = 0.68 GeV^ is 
found. This corresponds to a value for hraax of 0.5 GeV~^, 
and is what we will use in the fits of the next section. 



VI. SPECIFIC FITS 

In this section we provide examples of evolved fits, ob- 
tained by following the steps of Sec. |V] with specific fits 
for the input distributions. We remind the reader that 
our numerical calculations correspond to the Sivers func- 
tion of SIDIS, and that they acquire an overall minus sign 
in the Drell-Yan process. 



A. Bochum Fits 

The fits of Rcf. [T^] use a Gaussian to describe the 
HERMES measurements '36'] which were performed with 
an average of 2.41 GeV^. We refer to these as the 
Bochum fits. The function corresponding to f^^J [x) in 



Eq. ^ is 



_ up/down 



IT 



(x) 



± 



2Mi 



Bochum 

The fit parameters are 



Ax^-'il-xf. (53) 



A^O.n, & = 0.66. 

In the Bochum fits, the parameter corresponding to {k'^)^ 
in Eq. (1511) is assumed to be independent of fiavor and 



lies between 0.10 and 0.32 GeV^. We take 



Bochum 



= (4>0 



Bochum 



= 0.2GeV^ 



(54) 



which corresponds to the "best fit" scenario of Ref. jl4j . 

Samples of the result of using the Bochum fits in 
Eq. (|44p to evolve to different Q are shown in the up- 
per panel of Fig. [1] The curves are shown for Q ~ 
\/2.4, 5, 91.19 GeV since these are also the values already 
used to illustrate the evolution of the unpolarized distri- 
bution functions in Ref. [23 . 



B. Torino Fits 

Next we consider the fits of Ref. [ip which in- 
corporated data from both HERMES [sj and COM- 
PASS [H, [3^. Again, the scale for the initial distribu- 
tions is = 2.4 GeV^. We refer to these as the Torino 
jits. The function corresponding to f^^{x) in Eq. ([STjl is 



/iV(^) 

where 



Torino Mi(fcy) 



^^'^^'Uf{x)Mx){kl}o, 



Mf{x) ^ Nf x'^f (1 - x^f ^ '^7/.^ , 



(55) 



(56) 



and // (x) is the unpolarized parton distribution function 
for quarks of fiavor /. The fit parameters Nf, af,/3f are 



iVu = 0.35, 
iVd 0.90, 



Qfu — 0.73, 
ad = 1.08, 



I3u = 3.46, (57) 
/3d = 3.46, (58) 



and = 0.34 GeV^, (k^) 0.25 GeV^. The Gaussian 
slope parameter of the initial input distribution in the 
Torino fits is again flavor-independent and is 



t)o Torino 



= (4)o 



Torino 



Miikl) 



Ml 



{kl) 



(59) 



For the integrated PDFs in Eq. we have used the 

lowest-order MSTW parametrizations [40l - l43j . Samples 
of the evolved Torino fits are shown in the lower panel of 
Fig.ffl 

Note that there is over a factor of 2 difference between 
the Torino and the Bochum fits, and this gives a rough 
indication of the uncertainty involved in current treat- 
ments. A discussion of the difference in the two methods 
can be found in Ref. (i^ . 

We do hope for future improvements of the fits. A 
very recent parametrization of the nonperturbative in- 
put was presented in Ref. [1^. The results are similar 
to the Torino fits above, but utilize a relation to gener- 
alized parton distributions, and allow for a connection 
to a quantification of parton angular momentum. Mor- 
ever, model calculations, such as in Refs. [4a, |43, and 
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Up Quark Sivers Function 
x = 0.1 




(GeV) 

FIG. 1: (Color online.) The (negative of the) up quark Sivers function at x = 0.1 evolved from Q — V2.4 GeV(solid maroon) 
to Q = 5 GeV(dashed blue) and Q = 91.19 GeV(dot-dashed red). The upper plot is found by evolving; the Gaussian fits of 
the Bochum group 14] and the lower plot is found by evolving the Gaussian fits of the Torino group [ISl]. In the case of the 
Bochum fits, the down quark Sivers function is just the negative of the up quark one. For the Torino fits, the down quark 
Sivers function is obtained by multiplying the up quark Sivers function by —1.35. These functions acquire an overall reversal 
of sign if used in Drell-Yan. 



lattice QCD calculations [48| can aid in providing mean- 
ingful parametrizations of the nonperturbative input over 
the whole of phase space and open up interesting ques- 
tions regarding the matching of purely nonperturbative 
descriptions of the Sivers function to pQCD. 

C. Evolved Gaussian Parametrizations 

Figure [T] suggests that, apart from the tail at large 
/ct, the Sivers function continues to be well described by 
a Gaussian shape, even after evolution to large Q. To 
describe the evolution of a purely Gaussian parametriza- 
tion, with the x and fcr dependence factorized, requires 
only a specification of the scale dependence of the Gaus- 
sian parameters. This saves having to directly calculate 
Eq. (|44|) . and its transformation to momentum space, 
separately for each value of Q and x. Because of the 
general convenience of working with Gaussian functions, 
we have obtained Gaussian fits for a range of Q starting 
at Q = \/2A GeV for the Bochum and Torino fits up 
to Q = 90 GeV. The fits are obtained using the Wol- 
fram Mathematica 7 FindFit routine, and examples 
are shown as the dashed curves in Fig. [2] A table of the 
resulting values for the Gaussian parameters is shown in 
Table H] (Fortran, C-|— 1-, and Wolfram Mathematica 
7 code that produce evolved Gaussian fits is available 



at m.) 

In Fig. [21 we illustrate the quality of the Gaussian 
fits to the Sivers function at intermediate and large 
g (Q = 5 GeV and 91.19 GeV, respectively). In 
practice, the Sivers effect is often probed via observ- 
ables like Eq. ([52]) . so we have plotted the integrand, 
—2TTk^F^rp^^{x,kT;n,Q)- Note that, after the evolution 

to large Q, the —2TTkj,F^^^{x,kT;fJ,,Q) acquires a very 
broad tail for both the Bochum and Torino fits. The 
tail falls off slowly; for Q = 91.19 GeV, the ratio of the 
value of the Bochum fit at fcy = 10 GeV to the value at 
kx = 5 GeV is about 0.65. This is roughly consistent 
with the l/kx fall-off at large kr that is expected from 
the power counting arguments in Sec. lIIICl The last two 
columns in Table U show the values of kx where the ra- 
tio of the Gaussian fits to the original Sivers functions 
is 0.8. That is, above fc!^™™ (GeV) the Gaussian fits to 
the evolved Torino Sivers function drop to less than 0.8 
of the original evolved Sivers function and similarly for 

jLBochum 
"r.max ■ 

That the description at small kx remains Gaussian is 
not entirely surprising given that the input we use for 
the nonperturbative evolution is Gaussian (gKibr) (x b^). 
However, it should be emphasized that the perturbative 
contribution to evolution results in a substantial modifi- 
cation of the shape and normalization of the TMD PDF, 
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TABLE I: Table of evolved Gaussian parameters, obtained by fitting Gaussians to the evolved Bochum and Torino fixed-scale 
fits. The fits are for xA^ Fj:^pf(x,kT; ^j,,Cf) and are related to F^-^ {x, kr; ^, Cf) via Eq. H46|) . The parameters are listed 
for the up quark distributions at x = 0.1; the Sivers function at different values of x can be found by multiplying by the 
appropriate ratios obtained from Eqs. (|53II55|I . The Gaussian slope parameter is the same for the up and down quarks. The 
normalization parameters a^p are related to the down quark normalizations by fldown"™ = — oiup'^'^"™ and ffldown" ~ — LSSaJf™'"". 
The last two columns, k^°^^^ and fcy°maxi are the values of fcr above which the Gaussian fits drop to less than a ratio of 0.8 
of the Sivers functions calculated directly from Eq. (1441) . 
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even at low kx- Therefore, Table |T] is not the result of 
simply Fourier transforming the nonperturbative contri- 
bution to Eq. (gH). Rather, to get the right TMD PDF, 
even when using a Gaussian approximation for low kx, 
the full pQCD evolution must be included. We find that 
difference between the fitted Gaussian and the result ob- 
tained by naively Fourier transforming the nonperturba- 
tive part of the evolution is similar to what was found for 
the unpolarized TMD PDF (see Fig. 2 of Ref. [H). 

The presence of the tail illustrates the danger in eval- 
uating moment integrals like Eq. (j52p without a careful 
account of the large-Zcy behavior. For Q — 91.19 GeV, 
there is more than 40% suppression in the integral of the 
curves in Fig. [2] from to 10 GeV when the Gaussian fit 
is used rather than the fit including the tail. (Note that 
in principle the integral should be extended to order Q.) 
For the Q = 5 GeV curves, integrated up to 5 GeV, the 



corresponding suppression is only about 9%. 

By contrast, at low-Zc^ the Gaussian functions, shown 
as the dashed curves in Fig. [21 provide excellent approx- 
imations to the evolved Sivers function. This suggests 
that the evolved Gaussian approximation is especially 
suited to low-Q/low-ZcT studies. A sample of evolved 
Gaussian fits for lower Q is shown in Fig. |31 



VII. DISCUSSION AND CONCLUSIONS 

Many of the recent phenomenological efforts related 
to the study of transverse polarization effects in TMDs 
have assumed a lowest-order, generalized parton-model 
(GPM) picture (50| and work within a rather narrow 
range of energy scales. However, the full power of factor- 
ization theorems lies in their ability to make predictions 
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FIG. 2: (Color online.) The up quark Sivers function at Q = 5 GeV and Q — 91.19 GeV (solid curves) and the corresponding 
Gaussian fit for the low-kr region (dashed curves). Note that the function plotted is the Sivers function multiplied by — 27rfcy. 
The upper panel is obtained by evolving the Gaussian fits of the Bochum group [ij and lower panel is obtained by evolving the 
Gaussian fits of the Torino group Jjj] . Below each plot, the ratio between a Gaussian fit and the evolved function including 
the tail is also shown. 



for a variety of processes over a wide range of energy 
scales. In this article, we have explained the steps for 
implementing evolution for polarization dependent TMD 
PDFs, specifically illustrated with the Sivers function. 
The basic method is the CSS formalism [l^, [13, [3l| . 
with the specific formulation of the TMD-factorization 
formalism given recently in Ref. [2l|. An advantage of 
the most up-to-date TMD-factorization formula is that 
it is written in a form closely analogous to the GPM 
(see Eq. ([1])), with explicit definitions for the individ- 
ual TMDs. Therefore, existing treatments that rely on 
a GPM framework need only to replace the unevolved 
TMDs with the evolved ones. An important aspect of our 
approach is that it relies on a genuine, complete TMD- 



factorization formalism, to be contrasted with the resum- 
mation methodology that has often been relied on in the 
past to treat many aspects of TMD physics. That is, the 
TMD-factorization formalism provides, from the outset, 
a consistent treatment of factorization for the full range 
of kx (or, equivalently, the full range of bx in coordinate 
space). 

Fortunately, many of the results obtained from the 
treatment of unpolarized TMDs can be carried over di- 
rectly to the polarization dependent case, including the 
calculation of the anomalous dimensions 7f, 7d and "/k, 
and the CS evolution kernel K, in both its calculable per- 
turbative part and its nonperturbative part gxibx) that 
is known from fits to unpolarized Drell-Yan. An impor- 
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FIG. 3: (Color online.) The evolving Gaussian parameters for —2TTk^F^'^{x,kT;fJ.,Q) for a range of Q obtained from the 
Torino and Bochum fits. Table |I] lists the Gaussian parameters for a selection of Q. 



tant difference from the unpolarized case is in the match- 
ing at large-fc^. In the unpolarized case, the TMD PDF 
(or FF) matches to a twist-2 coUinear factorization treat- 
ment at large fc^, whereas the Sivers function matches 
to a tvifist-3 collinear factorization treatment related to 
the Qiu-Sterman formalism, as in Eq. (j47l) . Thus, the 
treatment provided in this article unifies several different 
aspects of TMD physics. 

It is vsrorth commenting on the often repeated state- 
ment (see, e.g., Ref. [Slj) that calculations in covariant 
gauges are impractical or inconvenient, and that working 
in light-cone gauge is therefore preferred. In our work, 
we find that the opposite is true. Namely, the calculation 
of the perturbative parts (at least to order as) follows 
clear-cut steps in Feynman gauge, while the derivation 
of TMD-factorization theorems is much more direct in 
Feynman gauge than in light-cone gauge. (Indeed, we 
are not aware of the existence of a detailed light-cone 
gauge derivation of TMD factorization.) Moreover, once 
the calculation of the perturbative parts has been per- 
formed in Feynman gauge, a generalized parton-model in- 
terpretation follows directly from the TMD-factorization 
formula in Eq. ([T}. For these reasons, we advocate con- 
tinuing to work in Feynman gauge for both calculations 
and derivations. 

We have implemented the evolution explicitly using 
as input the already known 7^?, and (supplied 
for easy reference in the Appendix, previous fixed-scale 
Gaussian fits of the Sivers function at low-Q 14, j^, and 
previous fits of the CSS formalism to DY |33| . For the ex- 
plicit calculations in the present article, we have focused 
only on the low-fcr region where we need not be con- 
cerned with the treatment of the Qiu-Sterman formalism 
at large fc^, and the approximations of Sec. IVl make sense. 



The resulting evolved momentum-space Sivers functions 
are shown in Fig. [T] Comparing with Fig. 1 of Ref. [l^l 
for the evolution of the unpolarized TMD PDF, one sees 
even more suppression as Q is increased than in the un- 
polarized case. Also note that a significant perturbative 
tail is generated at large Q as shown in Fig. [21 We reem- 
phasize that this should be kept in mind when evaluating 
integrals like Eq. ((52|) . 

Gaussian parametrizations are particularly convenient 
for doing explicit calculations. Therefore, we have tested 
the quality of Gaussian fits after evolution to large Q 
and find that the Gaussian function provides an excellent 
approximation to the Sivers function at small fcr, even 
for Q ss 90.0 GeV. We have made these fits available, as 
well as code for generating evolved TMDs at a website 
maintained by two of us (Aybat and Rogers) (49| . 

Much work remains to be done in the effort to connect 
a full QCD treatment of TMDs with phenomenology. An 
explicit implementation of the matching to the twist-3 
Qiu-Sterman formalism is still needed, and will be partic- 
ularly important for a correct treatment of fcT-weighted 
observables in which the extra fc^ factors enhance the 
contribution from the large fc^ region. The recent work 
of Ref. [m may help. Moreover, as new data become 
available for both polarized and unpolarized cross sec- 
tions, it will be useful to construct new fits that include 
evolution from the beginning. Finally, explicit calcula- 
tions, analogous to the ones presented here, need to be 
applied to the other TMDs like the Boer-Mulders and 
Collins functions. 

At large Q, the shape of the distribution is especially 
sensitive to the value of 6max, 52 and the functional form 
of (6t)- Reference [s^, for example, finds that a larger 
value of bmax is preferred, along with a corresponding 
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change in 52 • Furthermore, Refs. [12, find advantages 
to using a different functionai form, ~ iP'/^ rather than 
~ for gKibr)- This should be taken into account in 
future improvements to the fits. The particular set of 
parameters used in the calculations in the present article 
were chosen both because of their simplicity and because 
they correspond to the current state-of-the-art of global 
fits to the unpolarized Drell-Yan cross section. 

In the future, model calculations (see, e.g., and ref- 
erences therein for an overview) can be potentially help- 
ful for fixing nonperturbative input. Certain models also 
lead to nonperturbative input distributions that deviate 
from the Gaussian ansatz. Conversely, incorporating evo- 
lution into model calculations can help establish the scale 
appropriate to the model. 

Theoretical uncertainties in the TMD fits, both for un- 
polarized and polarized TMDs, can be reduced by in- 
cluding higher-order results for the anomalous dimen- 
sions and the CSS kernel K (in the perturbative region). 
Fortunately, as we have discussed in this paper, these 
anomalous dimensions and the kernel K are the same for 
unpolarized TMDs and the Sivers function. Therefore 
by calculating them at next-to-next-to-leading order in 
pQCD, we can reduce the theoretical uncertainties for 
both unpolarized and polarized TMDs at the same time. 

The ultimate goal is to obtain sets of TMD PDFs and 
FFs that can be used in a way that is closely analogous 
to what already exists for processes that use collinear 
factorization. Namely, we would like to obtain a set of 
TMD fits based on precise TMD definitions such that 
they can be reliably used to make predictions. 
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APPENDIX: ANOMALOUS DIMENSIONS ETC. 

Here we list the MS-scheme anomalous dimensions [2l| 
that were used in, for example, Eqs. and (|T7)) : 

7f(/^; (r/f^') = a,^ Q - In (^^)) + 0{al). (60) 

At order a^, the quark TMD FF anomalous dimension 
is the same as for the TMD PDF. The CS kernel up to 
order as in space is 

Ki^i, hr) = [IiV&t) - ln4 + 27e] + 0{al). 

n 

(61) 

The anomalous dimension of K is up to order as, 

7K(M)=2^+0(aD- (62) 
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